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Abstract 

During the last ice age there were several quasi-periodic abrupt warming events. The 
climatic effects of the so-called Dansgaard-Oeschger (DO) events were felt globally, al¬ 
though the North Atlantic experienced the largest and most abrupt temperature anoma¬ 
lies. Similar but weaker oscillations also took place during the interglacial period. This 
paper proposes an auto-oscillatory mechanism between sea ice and convective deep water 
formation in the north Atlantic as the source of the persistent cycles. A simple dynamical 
model is constructed by coupling and slightly modifying two existing models of ocean cir¬ 
culation and sea ice. The model exhibits mixed mode oscillations, consisting of decadal 
scale small amplitude oscillations, and a large amplitude relaxation fluctuation. The 
decadal oscillations occur due to the insulating effect of sea ice and leads to periodic 
ventilation of heat from the polar ocean. Gradually an instability builds up in the polar 
column and results in an abrupt initiation of convection and polar warming. The un¬ 
stable convective state relaxes back to the small amplitude oscillations from where the 
process repeats in a self-sustained manner. Freshwater pulses mimicking Heinrich events 
cause the oscillations to be grouped into packets of progressively weaker fluctuations, 
as also observed in the proxy records. Modulation of this stable oscillation mechanism 
by freshwater and insolation variations could account for the distribution and pacing of 
DO and Bond events. Physical aspects of the system such as sea ice extent and oceanic 
advective flow rates could determine the characteristic 1,500 year timescale of DO events. 

Keywords: Dansgaard-Oeschger events; Heinrich events; sea ice; deep water formation; 
climate oscillations; box model 


1. Introduction 

Climate proxy data from ice (?) and sediment (?) cores reveal that the climate of the 
north Atlantic underwent large quasi periodic fluctuations during the last glacial period 
(figure]^. These are the so-called Dansgaard-Oeschger (DO) events which occurred 
in intervals of about 1,500 years. Fluctuations on a similar timescale but with lesser 
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intensity are also seen during the Holocene warm period, known as Bond events (?). 
This paper proposes an auto-oscillatory mechanism between sea ice and convection in 
the north Atlantic to be the driving source of the pervasive millennial-scale cycles. A 
coupled sea ice-ocean box model is used to show the existence of stable oscillations whose 
pacing can be modulated by glacial meltwater pulses to produce temporal patterns very 
similar to the paleo-proxy record. 



Figure 1: Oxygen-18 isotopic concentrations in Greenland ice cores (NGRIP) showing the variation of 
near surface air temperatures over the last 100,000 years. Orange arrows mark individual DO events and 
the blue vertical stripes mark Heinrich events numbered 1 through 6, and the Younger Dryas flooding 
event. 


The DO events were characterized by an abrupt warming on a decadal scale followed 
by slower cooling. These abrupt warming events were felt most strongly in the north 
Atlantic. Proxy records from Greenland (?) and Antarctica (?) show that the tempera¬ 
ture anomalies in the northern and southern hemispheres were out of phase by about 800 
to 1,200 years, the magnitude of anomalies being much smaller over Antarctica. This 
‘bipolar seesaw’ behavior suggests a global reorganization of heat flow on the surface 
via changes in oceanic circulation (?). If the total amount of heat in the climate sys¬ 
tem is conserved then warming in one region must be balanced by cooling elsewhere on 
timescales comparable to the transit time of the global oceanic circulation (?). 

The warming over Greenland was also more abrupt than elsewhere, suggesting that 
the source of the heat anomaly was located in the north Atlantic. At present, a region 
of convective deep water formation exists in the Norwegian seas. The cooling and sub¬ 
sequent sinking of surface waters through buoyant instability provides a mechanism for 
the deep ocean to ventilate heat to the atmosphere. Freshwater perturbations via sudden 
bursts of glacial meltwater, i.e. Heinrich events, could disrupt the convective process and 
abruptly alter the overturning circulation and climate state (??). 

Heinrich events preceded some, but not all DO events. While there is little doubt 
that a large freshwater perturbation can disrupt and destabilize the convective process, 
it is not clear how a single pulse of freshwater could trigger multiple DO cycles. There is 
no evidence either to suggest that each DO event was triggered by a separate meltwater 
pulse. The timescale for ice sheet instability in models (?) is also signihcantly greater 
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than the pacing of DO events. Thus Heinrich events cannot be the driver of DO events. 

A number of hypotheses have been proposed to explain the source of oscillations. 
Some propose astronomical forcing through variations in solar output (?) or tidal res¬ 
onance from the moon’s orbit (?). Other studies suggest internal climatic mechanisms 
whereby steady salinity forcing, salting at low latitudes and freshening at high, can pro¬ 
duce self-sustained oscillations (??????). These mechanisms do not explain many of the 
observed characteristics of DO events. For example, proxies of solar output do not show 
correlation with the climatological anomalies from DO events (?). The lunar hypothe¬ 
sis, which connects tidal resonances to the destabilization of ice sheets and consequent 
disruption of the overturning circulation, do not explain the variability in the pacing of 
DO events through time. Internal thermohaline oscillations are unstable to freshwater 
perturbations and are unlikely to have persisted through the glacial and Holocene epochs 
and through large variations in surface freshwater conditions. A stable oscillation mech¬ 
anism is needed that is robust against freshwater anomalies and glacial-interglacial scale 
changes in the background climate. 

In this paper a sea ice-convection oscillator is proposed as the pacemaker for DO 
events. Some previous studies have noted the importance of sea ice in initiating polar 
convection (??). A Van der Pol type oscillator was proposed by Saltzman (?) in which 
sea ice cover regulates the positive feedback from oceanic uptake of CO 2 to produce 
relaxation oscillations. Although atmospheric CO 2 variations were unlikely to be the 
cause for DO events, as the convective region is not large enough to account for the 
required CO 2 flux, sea ice could in a similar manner regulate the heat content of the 
deep ocean (?). This insulating property of sea ice, its on/off presence, and geophysical 
constraints on its spatial extent could make it a suitable component of a stable oscillation 
mechanism. The following sections describe a simple dynamical model composed of sea 
ice and an ocean box model which interact to produce stable millennial scale oscillations. 
The simplified modeling approach makes it possible to capture the minimum number of 
processes that are necessary to produce oscillations and to analyze the model’s dynamical 
behavior. 


2. A simple dynamical model 

The model consists of an ocean box model (?) coupled to a thermodynamic sea 
ice model (?). The two components are among of the simplest available models of the 
respective processes. In the ocean model there are three meridional boxes arranged in 
two vertical layers. Each box has a uniform temperature and salinity which can change 
dynamically. Advective transport between adjacent boxes is driven by the horizontal 
pressure in the top layer and the net flow in the model is conserved. Sea ice forms or 
melts on the top polar box if the temperature goes below or above the melting point. 
The dominating influence of sea ice on the ocean circulation is through the regulation 
of heat exchange between the top polar box and the atmosphere. There is no dynamic 
atmosphere, just a static temperature gradient applied to the top layer boxes. 

The following modifications are made to the original models: 

• The masses of individual ocean boxes are computed run-time and not kept hxed 
as in (?). This is done to avoid an unphysical density anomaly arising out of the 
convective mechanism. The mixing of two equal density water masses produces a 
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slightly higher density of the combined water mass and this one time-step anomaly 
destabilizes the system to produce the reported ‘free oscillations’. This scenario is 
an artifact of the discretized description of the model and not due to a physical 
process. With the mass correction, the sea ice decoupled ocean model does not 
exhibit oscillations. While there is more than one way to resolve this computational 
issue, such as by increasing the convective mixing time-step or by using higher order 
non-linear terms in the equation of state, updating the mass of the boxes run-time 
is the least invasive of these ways in preserving the fundamental aspects of the 
original models. 

• Diffusive mixing is not considered as its inclusion only affects the response timescale 
and not the qualitative character of the model’s behavior. 

• The extent of the southernmost box meridional zone is shortened to 45°- 60° N to 
only consider the northern arm of the overturning circulation. 

• The distribution fractions for precipitation is different owing to the altered geom¬ 
etry of the system. 

• Sea ice is assumed to be a perfect insulator, unlike (?) where it had non zero heat 
permeability. 

• Heat and salinity adjustments due to brine rejection and latent heat of fusion of ice 
is ignored. The contribution of these terms in the governing equations is effectively 
zero in one cycle of sea ice advance and retreat and their exclusion does not alter 
the qualitative character of the model. These terms are left out so that only the 
ingredients necessary for oscillations can be examined. 

A description of the basic set up of the model is given below which includes the 
modihcations to the original components. The default parameter values and units are 
listed in Table [2 

2.1. Geometry 

The ocean boxes are arranged in three meridional zones and two vertical layers. The 
same indexing scheme is used as in (?). Areas of the boxes are computed by treating the 
entire ocean as part of a shell enveloping a sphere of the Earth’s radius. Flows take place 
between boxes sharing a common face. The extent of sea ice is limited to the top polar 
box only. A schematic of the model geometry with the box indexing scheme is shown in 
hgure 12 

2.2. Governing equations 

The system is dehned by a set of ordinary differential equations accounting for the 
rate of change of heat and salt content of each box, and the rate of change of sea ice 
extent. 



Pf + Po4’i,jSj + Cf 


( 1 ) 


dS, 


( 2 ) 


df 


dt 



T 
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( 3 ) 
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Acceleration due to gravity 

9.8 m s ^ 

Cp 

Specific heat capacity of water 

4000 J kg-i K-i 

Po 

Reference density of sea water 

1028 kg m-3 

Pice 

Density of ice 

920 kg m“^ 

a 

Thermal expansion coefficient 

1.7x10"'* K"* 

P 

Haline contraction coefficient 

0.76 

To 

Freezing point of water 

273 K 

Lf 

Latent heat fusion of ice 

3.34x10'^ J kg-* 

A 

Heat exchange coefficient between atmosphere and ocean 

30 W m-2 K-* 

Tr 

Reference temperature 

283 K 

Sr 

Reference salinity 

0.035 

At 

Integration time step 

0.25 years 

V 

Thermal gradient 

-0.74 K/° lat 

e 

Net evaporation from surface tropical box 

12 mm/day 

qt 

Salinity forcing ratios for surface boxes 

-fO.3, -0.2, -0.1 

Tf 

Atmospheric temperature above box 1 

297 K 

c 

Advective transport coefficient (inverse of friction coefficient) 

10"® m^ s kg"* 

dice 

Thickness of sea ice 

2 m 

T 

growth-melt timescale of sea ice 

3 years 


Table 1: Description of model parameters. 


Here i is the index of a box. Terms on the right hand-side of equations Q and ([^ 
are contributions to heat and salinity fluxes via atmospheric forcing {F^ ,Ff), advective 
flows and convective mixing ^f). The flow terms, represent flows between 

boxes i and j, and i^ijTj and are implied summations over j to account for the 

net flux of heat and salinity into box i from its adjoining boxes. The convective mixing 
terms are temperature and salinity adjustments on uniformly mixing the top and bottom 
boxes in a column. 

In equation ^ f represents the fractional sea ice extent, ki is the fractional growth- 
melt rate constant of sea ice per unit difference in temperature from freezing, r is the 
timescale for sea ice growth-melt, and k 2 accounts for precipitation on existing sea ice. 
Sea ice also cuts off heat and salinity exchange between the top polar box and the 
atmosphere in proportion to /. 

2.3. Thermal forcing 

Each of the top layer boxes are made to relax to a static atmospheric temperature, 
with the air temperatures decreasing linearly with latitude. The heat exchange rate for 
boxes 1 and 2 are 

Ff = \{Tf - T,)A, (4) 

where Ai is the surface area of the box in contact with the atmosphere, Tf is the atmo¬ 
spheric temperature, and A is the heat exchange coefficient. For the polar box the area 
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Figure 2: Schematic of the box model. The solid and dashed arrows indicate up and downstream flows 
respectively. The height and width of the boxes are shown to scale, however the volumes of boxes is not 
to scale with the areas as shown due to the three dimensional geometry of the system. 


of contact is adjusted according to the fraction of area covered by sea ice 

F3^ = A(T“-T3)A3(1-/) (5) 

The atmospheric temperatures are prescribed by the thermal gradient parameter 77 , 
and a prescribed T“. 

n = v{e,-e,) + T^ ( 6 ) 

where 9i is the middle latitude of box i. 

2.4- Salinity (freshwater) forcing 

In addition to thermal forcing from the atmosphere, each top layer box receives 
salinity forcing due to evaporation or precipitation. Salinity flux from the atmosphere is 
positive for net evaporation, as in box 1, and negative for net precipitation in the other 
two boxes. The addition of salinity is such that the total amount of salt is conserved in 
the system. The distribution fractions of precipitation is different from (?) to scale for 
evaporation and precipitation in the absence of a tropical ocean box. The distribution 
fractions are 


qi — +0.3 

( 7 ) 

92 = -0.2 

(8) 

93 = - 0.1 

( 9 ) 

The salinity flux terms for each top layer box is given by 

Ff = poSoeAiqi,i = {1,2} 

(10) 

F 3 = PoSoeAi{l - f)q3 

(11) 


where po and So are reference density and salinity, e is the net evaporation from box 1 
in mm/day. 
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2.5. Advection 

Advective flows in the model are driven by horizontal pressure gradients given by 




( 12 ) 


where the parameter C is the inverse of a friction coefficient whose value chosen to 
give appropriate flow rates for the North Atlantic, is the area of contact between 

adjoining boxes, and are the respective pressures in the boxes. Pressure is 

formulated in such a way that the net pressure integrated along a vertical column is zero. 


Pi — „ 7 I 7 {hiPi + hi_|_3Pi_|_3)g 

z rii -i- 

hi 

Pi+3 = -W^^2 

^2 + 3 


(13) 

(14) 


where hi and /ii +3 are the depths, and pi and pi +3 are the densities of the top and bottom 
boxes in a vertical column. Density is computed from the linear equation of state given 
as 

p,=Po{l-a{T,-Tr)+^{S,-Sr)) (15) 

where a and j3 are thermal expansion and haline contraction coefficients respectively, 
and Tf, are reference temperature and salinity values. 

Volume conservation requires flows between any pair of top layer boxes to be matched 
by flow in the opposite direction in the bottom layer. The up/downwelling flux between 
boxes 2 and 5 is similarly derived from volume constraints as 


02 « = 05,2 = 02-01 


(16) 


A positive 0^ implies poleward flow on the surface, and a positive 0^, implies up- 
welling. 

2.6. Convection 

The convective scheme mixes the heat and salt contents uniformly between the boxes 
in a vertical column when a buoyant instability, pi > pi+ 3 , occurs. The equalized tem¬ 
perature and salinity values are given by 


5,; = 


T, = 


T ^^2+3 *^2+3 
m* -I- mi+3 

miTi + 

m* -I- mj+3 


So that the heat and salinity adjustment terms in equations Q and Q are then 


= 


= 


T^-T, 
_ At 
S^-S, 
At 


(17) 

(18) 

(19) 

( 20 ) 


where At is the time step over which convection is carried out, which by default is equal 
to the integration time step. 
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2.1. Sea ice 


Sea ice forms/melts when the temperature of the surface polar box goes below/above 
freezing. The governing growth rate by volume (?) is given by 


^bice 

dt 


^^^%^(To-r3) + /p3 

AceT/T 


( 21 ) 


where V 3 is the volume of the polar box, r is a prescribed timescale for sea ice growth, 
Tq is the freezing temperature, and is the precipitation rate by volume per unit time 
given by 

P3 = A^elqsl (22) 

so that any precipitation on existing sea ice is turned into additional sea ice. The volume 
of ice formed is spread over the polar box with thickness dice- The fractional growth rate 
in sea ice is then 


# ^ 1 dVice 

dt dicers dt 

^ice-^3 y/^ice^/'^ 

which is expressed in a more compact form as equation ([^. 
then given by 



(23) 


The constants ki and ^2 are 


ki 

k2 


1 PoCpVs 
dicers Piced^f 
P3 

dicers 


(24) 

(25) 


As discussed earlier the only influence of sea ice in this model is the alteration of 
heat and salinity exchange of the polar box with the atmosphere. Brine rejection and 
latent heat of fusion are ignored as their contributions in one cycle of sea ice advance 
and retreat equal zero. Although those terms affect the timescale of oscillations, their 
exclusion does not change the qualitative character of the model and they are therefore 
ignored in order to focus on the driving ingredients only. 


3. Model dynamics 

The model exhibits two distinct equilibrium states and two self-sustained oscillating 
states. The steady states have opposite flow directions, termed as thermal (TH) for 
poleward flows in the top layer and haline (HA) in the opposite direction. One of the 
oscillating states (OS-2) is composed of two alternating modes, a small amplitude multi- 
decadal scale oscillation and a large millennial scale fluctuation (figure . The other 
oscillating state (OS-1) is the small amplitude oscillation mode sustained indefinitely. 

The large amplitude component in oscillations is initiated when the density of the 
top layer polar box exceeds the density of the bottom box. Convection occurs when the 
instability condition is met and it results in an abrupt increase in poleward advective 
flux and with it an abrupt polar warming. The convective state is transient following 
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Figure 3: Variation in advective fluxes i/ii (blue), r/)2 (purple), and sea ice (orange) with time. Sea ice 
extent is typically either on or off. 


which the system relaxes back to the small amplitude oscillations. During this phase 
sea ice and the polar advective flux, ' 02 , oscillate antagonistically with the advective flux 
switching between positive and negative, figure (|^. The amplitude of these oscillations 
increase with time until a convective instability is reached in the polar column. 

In each phase of sea ice advance and retreat the top polar box loses net heat. This is in 
spite of the insulating effect of sea ice. While sea ice covers the polar box its temperature 
increases until this causes the ice to melt and retreat. The lowering of density due to the 
temperature increase drives the advective flux poleward, and poleward heat transport. 
The high temperatures also mean higher ocean-atmosphere temperature gradient, which 
results in a greater loss of heat. 

High mixing rates in the convective phase remove the density anomalies that resulted 
in the convective instability. This allows the system to relax back to the small amplitude 
oscillations, which in the absence of sea ice would be a stable fixed point. The existence of 
sea ice destabilizes the HA equilibrium. For appropriate strengths of the salinity forcing, 
the small oscillations are interrupted when the system touches the convective manifold 
and gets forced away from the HA fixed point. A three-dimensional phase space that 
best captures the dynamical behavior of the system is spanned by the top and bottom 
polar densities and sea ice extent. Figure (§ shows one entire oscillation cycle, consisting 
of both the small and large amplitude components. 

The model’s response is systematically examined in a large range of thermal and 
freshwater forcing values. Independent simulations with identical initial conditions are 
carried out in 50x50 pairs of forcing values. Detection of oscillations and/or steady states 
from each of the generated time series is automated. Equilibrium flux values, oscillation 
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Figure 4: Antagonistic variations in sea ice extent and These decadal scale oscillations are precursors 
to a convective instability in the polar water column. 



Figure 5: A three dimensional phase space spanned by the densities of the polar boxes and sea ice extent. 
When the small amplitude high frequency oscillations touch the convective manifold, where p 3 = pg, an 
abrupt transition to a thermal mode occurs. 
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frequencies and amplitudes are recorded for each simulation. Figure (§ shows the states 
in the forcing space. For the set of simulations where sea ice is decoupled from the ocean 
model, the state space only has TH and HA states and no oscillations, hgure Q. 
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Figure 6: Model states in the two-dimensional forcing space. Larger magnitudes of r] imply steeper 
atmospheric temperature gradients. 


In both the sea ice coupled and decoupled cases, TH states occur when the relative 
strength of the thermal forcing is greater than the salinity forcing. HA states are only 
observed for extremely large salinity forcing where sea ice is coupled. OS-2 states in the 
sea ice coupled model occurs along the edge of TH and HA states where the relative 
strengths of the two forcings nearly balance each other. 

3.1. Bifurcations 

The bifurcation structure is examined by slowly varying the salinity forcing parameter 
and observing the model’s behavior, while keeping the thermal forcing parameter fixed. 
For every new parameter value the model is allowed to run for a sufficient length of time 
to let it come to equilibrium or sustained oscillations, as the case may be. Traversing 
through ascending values of e shows all the four modes, TH, OS-2, OS-1, and HA in 
successive order, figure Q. 

A steady TH state near the bifurcation point to OS-2 is given a series of pertur¬ 
bations by adding freshwater to the top polar box. Small magnitudes of perturbation 
decay exponentially and the system spirals and settles back into the steady state. If the 
perturbation is larger than a certain threshold, figure the trajectory diverges far, 
into the haline regime. This behavior is characteristic of subcritical Hopf bifurcations. 
As expected, the system exhibits hysteresis if the direction of varying e is reversed. 

The OS-2 state is a mixed mode, with two oscillation components. The small ampli¬ 
tude high frequency component is near an unstable HA fixed point. It is the coupling of 
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Figure 7: Only steady TH and HA states occur in the ocean model decoupled from sea ice. 



Figure 8: Bifurcation diagram for slowly varying salinity forcing parameter e. 
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Figure 9: Subcritical Hopf bifurcation from TH to OS-2 state. A series of progressively larger pertur¬ 
bations are applied to the steady TH state near the bifurcation point. Perturbations with amplitudes 
below a threshold value decay exponentially in time, while those larger are repelled to the haline regime 
of flow. 


sea ice that makes this fixed point unstable, as without it the small oscillations do not 
occur. However, before the oscillations can settle on a limit cycle, they could encounter 
the convective manifold and get repelled towards a transient thermal mode. The abrupt 
jump and relaxation back to the haline mode is the large amplitude low frequency com¬ 
ponent of OS-2. Higher e values push the limit cycles away from the convective manifold, 
so that the convective detour never takes place and the system persists in the small am¬ 
plitude limit cycle mode, or OS-1. Thus the bifurcation from OS-2 to OS-1 is of the 
homoclinic type. 

The OS-1 to HA bifurcation occurs when sea ice completely fails to grow. As soon as 
sea can grow a small amplitude limit cycle appears, which suggests that a supercritical 
Hopf bifurcation has occurred. 

The dynamical properties of the ocean model by itself is similar to that of Stommel’s 
box model in that it exhibits stable haline equilibrium modes and spirally stable thermal 
mode. ? prove that Stommel-like models cannot exhibit oscillations under steady forcing. 

4. Period dependence on parameters 

Oscillation timescales are determined by four main parameters - the ratio of the 
thermal to salinity forcing ry/e, the growth-melt timescale for sea ice r, the overturning 
flow rate constant C, and the maximum extent of sea ice. With the exception of the 
growth timescale, all other parameters varied significantly over the course of the last 
glacial period and the response of the sea ice-convection oscillator to those variations 
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may have played a significant role in determining the climate of the north Atlantic. The 
model’s response is analyzed for different values of these parameters. 


4-.1. Ratio of thermal to salinity forcing magnitudes 

The mixed mode oscillations occur when the effects of thermal and salinity forcing 
nearly balance each other. Periods of oscillations have two components, the timescale 
for a convective instability to build up in the polar column, and the relaxation timescale 
for the system to return from the transient thermal mode. Strong atmospheric thermal 
gradients increase the attractive basin for thermal circulation modes and therefore tend to 
increase the relaxation timescale, which has the most pronounced effect on the oscillation 
period. Figure (101 tabulates the periods of all OS-2 states from figure (1^. It shows 


that the periods vary from a few hundred to a few thousand years. 
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Figure 10: Distribution of periods over the ratio of forcing strength parameters. Larger negative values 
of the ratio imply greater relative strength of thermal forcing. 


While 1,500 years is the most commonly cited period for DO events, the proxy records 
show that their pacing varied quite significantly during the last 100,000 years. The most 
regular periodic occurrence of DO events is in the window of 30,000 and 50,000 years 
before present, which was also the time when high latitude summer insolation had the 
least amount of variation. In other time windows the pacing of DO events appear to be 
modulated by insolation variations. The distribution and duration of abrupt warming 
events during the last glacial period was probably shaped in large part by a combination 
of solar and glacial meltwater forcing. 

4-2. Timescale for sea ice growth/decay 

The growth and melt timescale for sea ice depends on a large number of physical 
conditions in the real world, such as the spatial and seasonal variations in near surface 
temperatures, ocean-atmosphere heat fluxes and strengths of ocean currents. In the 
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model this timescale is prescribed to be in the 2-3 year range, which is largely an educated 
guess of how fast sea ice could grow over the Nordic seas under a static background climate 
similar to the last glacial period. During last glacial maximum the southern extent of 
winter sea ice was Iceland, whereas today it is near Svalbard. The area of sea between 
Svalbard and Iceland is considerably larger than the present seasonal variations in sea ice 
extent. Considering the annual seasonal cycle of growth and melt the default timescale 
for sea ice is reasonable. 

With all other parameters fixed at their default values the model is run for a range of 
sea ice timescales. An abrupt transition occurs around a value of 1.5 years, above which 
oscillations appear and approach periods of about 1,400 years (figure [TI|) . The threshold 
occurs where the insulating effect of sea ice in trapping heat in the ocean is matched 
by the dissipative heat loss by conduction and advection. The sensitivity for values of r 
above this threshold is low as only the high frequency component is affected and not the 
longer timescale relaxation mechanism. 
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Figure 11: Distribution of periods over the varying timescales of sea ice growth. No oscillations occur if 
r < 2 years. Above this value the sensitivity of periods to the growth timescale is low. 



4--3. Advective flux rate constant 

The parameter C prescribes the meridional flow rate per unit pressure gradient and 
its default value is based on current flow rates in the North Atlantic (?). However, this 
rate was likely very different between glacial and interglacial times (?). In the context of 
a relaxation oscillator, especially one where the flow is conserved, the rate of flow would 
influence timescales of both instability build up and relaxation. For small flow rates the 
heat build up in the ocean depths would fail to bring about convective instability, whereas 
high flow rates would equalize density anomalies efficiently and prevent oscillations. As a 
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result the model should exhibit oscillations only within a range of flow rates, as observed 
in simulations over a range of flow rates (figure [T^ . 
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Figure 12: Distribution of periods over the strength of advective flow rate. 


4 5 6 

Advective heat transport coefficient, 
C (x10“® nf s kg"^) 


4-.4- Maximum extent of sea ice 

The southern extent of sea ice in the north Atlantic is sensitive the background 
climatic state, glacial or interglacial, as evidenced by proxy data measuring sea ice extent 
(?). However among the other parameters that control the periodicity of the system, the 
sea ice extent is most constrained due to the positioning of land masses and structure of 
ocean currents in the north Atlantic. Could it be possible that the characteristic 1,500 
year period of DO events is due to the size of the ocean basin effected by sea ice? 

To test this hypothesis the area of the top polar box is varied between factors of 0.5 
and 5 of its default area, and oscillation periods observed in this range. Periods show high 
sensitivity to the polar area, which is also a measure of the maximum extent of sea ice. 


A peak in the periods is observed around a factor of 3, figure (13). The north Atlantic 


oscillator system in the real world could also have a similar non-linear sensitivity to the 
extent of sea ice. This would make certain oscillation periods favorable within certain 
windows of sea ice variance, and very different periods in other windows, as may have 
been the case between glacial and interglacial climates. 


5. Modulation by freshwater forcing 


The addition of freshwater to the ocean surface due to catastrophic collapse of ice 
sheets has been speculated to be the trigger for DO events. Simulations with intermediate 
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Figure 13: Oscillation periods for different polar areas. The horizontal axis is the scaling factor of the 
default polar area. 



complexity models do show that freshwater additions in the North Atlantic can abruptly 
destabilize and enhance the meridional circulation (??). The higher complexity models, 
however, do not seem to capture an internal climatic oscillator and consequently the 
model responses do not exhibit periodic behavior in the absence of periodic freshwater 
forcing. Typically in these model studies every DO event needs to be triggered by a 
separate freshwater pulse, or an underlying freshwater periodicity imposed on the system. 
The climate record, however, suggests an intrinsic periodicity within the climate system 
that was also sensitive to freshwater pulses. 

The effect of freshwater on the circulation state alone is examined by adding incre¬ 
mental amounts of freshwater to box 1, with no sea ice coupling on the polar box. The 
freshwater forcing is varied slowly, allowing the model to come to equilibrium in each step. 
A hysteresis is observed when the variation in freshwater additions is reversed, showing 
a bi-stable region and a stable HA mode, figure (141. In the absence of sea ice there is no 
mechanism for the circulation to undergo relaxation oscillations and consequently those 
are not observed. 

The effect of Heinrich events is simulated by constructing a corresponding freshwater 
signal and forcing the sea ice coupled model with it. The pulses are spaced in intervals 
of 7,500 years and are represented by a Gaussian function. A smaller magnitude skewed 
Gaussian is added in between pulses to represent increasing basal melting from ice sheets. 

The otherwise periodic response of the model is modulated by the freshwater forcing 
for Heinrich events, figure (15). Each Heinrich event triggers a series of progressively 
weaker abrupt warmings, similar to what is observed in the climate record. The gradual 
addition of freshwater to the surface in between Heinrich events decreases the meridional 
pressure gradient and results in weaker transient thermal modes of circulation. The 
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Figure 14: Steady states of the ocean box model decoupled from sea ice for varying rates of freshwater 
additions to the surface. A hysteresis is observed when the direction of freshwater additions is reversed. 


freshwater pulse causes the poleward circulation to abruptly cease, which initiates sea 
ice growth, which in turn causes a resumption of convection and poleward circulation. 

Variation of high latitude summer insolation during the last glacial period was least 
in the window of 35,000 to 50,000 years before present. This is the same time window 
where DO events occurred with the regularity simulated with the model, possibly due 
to minimal influence of insolation forcing in that period. Speculating further it could be 
stated that DO events were most pronounced during the glacial period, in comparison to 
the Holocene, because the background circulation was in a weaker state. A combination 
of solar and surface freshwater forcing could explain the climatic fluctuations during the 
last 100,000 years, which is out of the scope of this paper and largely a diagnosis of the 
system. 

6. Discussion and conclusions 

A simple dynamical model is used to show that sea ice and convection in the north 
Atlantic could interact to produce millennial scale quasi-periodic fluctuations. In the 
idealized model described here, the fluctuations are stable periodic orbits under a large 
range of climatic forcing. The timescale of model oscillations and the temporal pattern of 
variability is very similar to that of DO events and Bond Cycles. Time varying freshwater 
anomalies due to ice sheet discharge and variations in obliquity can modulate the intrinsic 
periodic response of the sea ice-convection oscillator. This oscillator could be at the heart 
of the millennial scale climatic variability known to have occurred during the last 100,000 
years. 
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Figure 15: Forcing by Heinrich events (blue curve) can modulate the response of the sea ice-convective 
oscillator into patterns that are similar to the temperature fluctuations recorded in ice core proxies. 


The oscillation is due a periodic build up of convective instability in the polar ocean. 
Sea ice cover, when present, insulates the ocean from losing heat to the atmosphere. 
However, over repeated cycles of sea ice advance and retreat the upper mixed oceanic 
layer undergoes a net loss of heat due to the large heat ventilation during retreat phases. 
The gradual increase in mixed layer density leads to convective instability which initiates 
an abrupt increase in vertical mixing in the polar water column. The high mixing rate 
removes the density anomalies and allows the circulation to relax back to its preferred 
state of small amplitude oscillations in sea ice advance and retreat. From here on the 
net positive heat loss from the mixed layer helps build up the convective anomaly and 
repeat the process in a self-sustained manner. 

Dynamically, the oscillations are composed of two components operating on a decadal 
and millennial timescale respectively. The presence of sea ice is responsible for the small 
amplitude oscillations, which are not observed in the ocean model decoupled from sea 
ice. These oscillations approach a limit cycle, but may be interrupted by the trajectories 
running into the convective manifold which force the system away from the unstable HA 
fixed point towards a transient thermal mode. The onset of convection brings about a 
high mixing rate which gradually removes the density anomalies and allow the circulation 
to relax back to the haline mode small oscillations. The mixed mode character of the 
sea ice oscillator means that there are at least three principal dynamical variables that 
are necessary ingredients. These are likely the horizontal and vertical density gradients 
in the north Atlantic ocean, and the extent of sea ice. Could it be that the well observed 
Atlantic Multi-decadal Oscillations are the small amplitude component of a mixed mode 
oscillation of the sea ice-convection system? 
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The period of oscillations in the model are dependent on the climatic forcing as well 
as physical attributes of the system. Sea ice extent is geophysically and climatologically 
limited within the corridor of Greenland-Icelandic-Norwegian seas and this could set the 
intrinsic timescale of 1,500 years. Freshwater pulses mimicking Heinrich events modulate 
the periodic response into repeating groups of successively weaker and shorter fluctua¬ 
tions. This response is similar to the Greenland temperature proxy signal in the window 
of 35,000 to 50,000 years before present. In this time window the summer insolation at 
high latitudes had very low variability as compared to the rest of the glacial period and 
thus surface freshwater forcing may have been the dominating influence. 

A box model, in spite of its simplicity, is an effective tool in identifying the fundamen¬ 
tal mechanisms and ingredients behind the climatic phenomenon addressed in this paper. 
The persistent fluctuations occurring through both glacial and interglacial periods hint 
at a robust oscillation mechanism for which the sea ice-convection oscillator is a suitable 
candidate. New mathematical techniques to study non-smooth dynamical systems such 
as this one are needed for further exploration of the dynamics of this climate oscillator. 
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